/*copyright 2009 Kevin Daley
#copyright 2009 Kevin Daley
#This file is a part of Cerebra.
#Cerebra is free software: you can redistribute it and/or modify
 #   it under the terms of the GNU  General Public License as published by
  #  the Free Software Foundation, either version 3 of the License, or
   # (at your option) any later version.

   # This program is distributed in the hope that it will be useful,
   # but WITHOUT ANY WARRANTY; without even the implied warranty of
   # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   # GNU  General Public License for more details.

   # You should have received a copy of the GNU  General Public License
   # along with this program.  If not, see <http://www.gnu.org/licenses/>.

#Kevin Daley, Lead Developer.  still.horse@gmail.com.*/

#define BLOCK_SIZE 16
#define kk 2*3.1415927f
#define lmb 1
struct params{
	float d, dt;
};
 typedef struct {
                   int width;
                   int height;
                   int stride;
                   float* elements;
               } Matrix;
__global__ void node 
( float* n  /*the number of instances of each descriptor, in row-major order*/,  float* spectrum /*the unique numeric identifier of each descriptor*/,  float* J /*the precomputed spectral Julia set, in row-major order*/,  float* j/*the boundary condition*/, int copy_point[2], float dt, struct params etc, int N, int bands){   				
	int2 index=make_int2((threadIdx.x+blockIdx.x*blockDim.x)%8+copy_point[0], (threadIdx.y+blockIdx.x*blockDim.x+copy_point[1])%8);
	int i=int(blockIdx.x*blockDim.x+threadIdx.x)/64;
	__shared__ float sdata[1024+80]; //256 vectors each for: old instance array, changes made to that, total energy, julia set
	float k;
	float2 delta;
	for(int lp=0; lp<N; lp++){
		k=spectrum[i];
		sdata[(threadIdx.x+threadIdx.y*BLOCK_SIZE)]=n[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x];
		n[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x]=0;
		sdata[255+BLOCK_SIZE+(threadIdx.x+threadIdx.y*BLOCK_SIZE)]=J[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x];
		if (sdata[255+BLOCK_SIZE+(threadIdx.x+threadIdx.y*BLOCK_SIZE)]!=0 or index.x%4+index.y%4==0){
		sdata[767+48+(threadIdx.x+threadIdx.y*BLOCK_SIZE)]=n[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x]*spectrum[i];
		delta=make_float2(-k*lmb*cos(dt*k*kk)*sdata[(threadIdx.x+threadIdx.y*BLOCK_SIZE)]/etc.d, lmb*k*sin(dt*k*kk)*sdata[(threadIdx.x+threadIdx.y*BLOCK_SIZE)]/etc.d);//find the displacement operator
		__syncthreads();
		n[int(i*(BLOCK_SIZE*gridDim.x)*(BLOCK_SIZE*gridDim.y)+index.x+index.y*BLOCK_SIZE*gridDim.x+delta.x+delta.y*BLOCK_SIZE)]+=(j[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x])*(-1/lmb*pow(k*n[int(i*256+index.x+index.y*BLOCK_SIZE*gridDim.x+delta.x+delta.y*BLOCK_SIZE)]/k*n[i*256+index.x+index.y*BLOCK_SIZE*gridDim.x],2.0f)*dt/spectrum[i]+1/lmb*pow(k*n[int(i*256+index.x+index.y*BLOCK_SIZE*gridDim.x+delta.x+delta.y*BLOCK_SIZE)]/k*n[int(i*256+index.x+index.y*BLOCK_SIZE*gridDim.x)],1.0f)*dt*spectrum[i]);
		}
	}
	}

__global__ void  calc_Julia(float* Jm, float* spectrum, float d){
                       //real and imaginary part of the constant c, determinate shape of the Julia Set
    float newRe, newIm, oldRe, oldIm;   //real and imaginary parts of new and old
    int maxIterations = 3000; //after how much iterations the function should stop

    //pick some values for the constant c, this determines the shape of the Julia Set
	
    //loop through every pixel
    int x=int(blockIdx.x*blockDim.x+threadIdx.x);
    int m=int(x/8);
    float cRe=spectrum[int(x/8)];
    float cIm=0;
    x=x%8;
    int y=int(blockIdx.y*blockDim.y+threadIdx.y);
    y=y%8;
        //calculate the initial real and imaginary part of z, based on the pixel location and zoom and position values
        newRe = x/d;
        newIm = y/d;
        //i will represent the number of iterations
	int i;
        //start the iteration process
        for( i = 0; i < maxIterations; i++)
        {
            //remember value of previous iteration
            oldRe = newRe;
            oldIm = newIm;
            //the actual iteration, the real and imaginary part are calculated
            newRe = oldRe * oldRe - oldIm * oldIm + cRe;
            newIm = 2 * oldRe * oldIm + cIm;
            //if the point is outside the circle with radius 2: stop
            if((newRe * newRe + newIm * newIm) > 4) break;
        }
	Jm[x+y*(gridDim.x*BLOCK_SIZE)+m*(BLOCK_SIZE*gridDim.x)*(BLOCK_SIZE*gridDim.y)]=(i < maxIterations);
}
		
	



               

               // Get a matrix element
               __device__ float GetElement(const Matrix A, int row, int col)
               {
                   return A.elements[row * A.stride + col];
               }
               // Set a matrix element
               __device__ void SetElement(Matrix A, int row, int col,
                                            float value)
               {
                   A.elements[row * A.stride + col] = value;
               }
               // Get the BLOCK_SIZExBLOCK_SIZE sub-matrix Asub of A that is
               // located col sub-matrices to the right and row sub-matrices down
               // from the upper-left corner of A
               __device__ Matrix GetSubMatrix(Matrix A, int row, int col)
               {
                   Matrix Asub;
                   Asub.width      = BLOCK_SIZE;
                   Asub.height     = BLOCK_SIZE;
                   Asub.stride     = A.stride;
                   Asub.elements = &A.elements[A.stride * BLOCK_SIZE * row
                                                         + BLOCK_SIZE * col];
                   return Asub;
               }

__global__ void MatAnd(Matrix A, Matrix B, Matrix C)
   {
       int row = threadIdx.y;
       int col = threadIdx.x;
	SetElement(C, row, col, GetElement(A, row, col) == GetElement(B, row, col));
   }
